Geog533 Lab 10

Question 1

Use R package UScensus2010county to complete the following tasks: (20 pt.)

library(UScensus2010county)

Question 1(a)

Plot a map of New York counties using the plot function.

data(new_york.county10)
shp <- new_york.county10
df <- shp@data
df

plot(shp)

Question 1(b)

Plot a map of New York counties using the qtm function.

install.packages("tmap")

Error in install.packages : Updating loaded packages

library(tmap)
qtm(shp)

Question 1(c)

How many counties in New York State?

df

print("There are 62 counties in New York State")

[1] "There are 62 counties in New York State"

Question 1(d)

What’s the 3-digit fips code of Broome County?

print("36007")

[1] "36007"

Question 1(e)

Compute descriptive statistics of the population column (P0010001), including total, minimum, maximum, mean, median, and skewness.

nypop <- df$P0010001
t <- sum(nypop)
t

[1] 19378102

a <- summary(nypop)
a

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   4836   51244   91301  312550  231060 2504700 

library(moments)
b <- skewness(nypop)
b

[1] 2.579801

Question 1(f)

Create a histogram and a boxplot of the population.

a <- hist(nypop)

a

$breaks
[1]       0  500000 1000000 1500000 2000000 2500000 3000000

$counts
[1] 53  3  3  1  1  1

$density
[1] 1.709677e-06 9.677419e-08 9.677419e-08 3.225806e-08 3.225806e-08 3.225806e-08

$mids
[1]  250000  750000 1250000 1750000 2250000 2750000

$xname
[1] "nypop"

$equidist
[1] TRUE

attr(,"class")
[1] "histogram"

b <- boxplot(nypop)

b

$stats
       [,1]
[1,]   4836
[2,]  51125
[3,]  91301
[4,] 234878
[5,] 468730

$n
[1] 62

$conf
          [,1]
[1,]  54429.09
[2,] 128172.91

$out
[1]  744344 1585873 2504700 1339532  949113 2230722  919040 1493350 1385108

$group
[1] 1 1 1 1 1 1 1 1 1

$names
[1] "1"

Question 2

Use R package UScensus2010tract to complete the following tasks: (20 pt.)

library(UScensus2010)
if(!require(UScensus2010tract)) install.tract("osx")
library(UScensus2010tract)

Question 2(a)

Plot a map of New York census tracts using the plot function.

data(new_york.tract10)
shp <- new_york.tract10
df <- shp@data
df

Restarting R session...

Question 2(b)

Compute the total population based on census tracts.

x

[1] 19378102

Question 2(c)

Select all census tracts in Broome County and plot the map.

BroomeCounty <- county(fips="007",name="Broome",state="NY",level="tract")
plot(BroomeCounty)

Question 2(d)

What’s the total population of Broome County?

Broomepop <- sum(BroomeCounty$P0010001)
Broomepop

[1] 200600

Question 2(e)

Create a histogram and a boxplot of population based on census tracts of Broome County.

hist(BroomeCounty$P0010001)

boxplot(BroomeCounty$P0010001)

Question 2(f)

Select the first five columns of the shapefile of Broome County census tract; add a new population ratio column (= census tract population / county population); save the new shapefile to the hard drive.

pr <- data.frame(BroomeCounty[,1:5])
pr$ratio <- BroomeCounty$P0010001/200600
pr

Question 3

Use R packages raster and ncdf4 to complete the following tasks: (20 pt.)

library(raster)
library(ncdf4)

Question 3(a)

Load the multi-band raster dataset “NDVI.nc” into R.

setwd("data/gimms_ndvi/")

Error in setwd("data/gimms_ndvi/") : cannot change working directory

Question 3(b)

Get the basic information about the dataset, including the number of rows, columns, cells, and bands; spatial resolution, extent, bounding box, and projection.

ndvi.rb <- brick("NDVI.nc")
ndvi.rb

class       : RasterBrick 
dimensions  : 360, 720, 259200, 36  (nrow, ncol, ncell, nlayers)
resolution  : 0.5, 0.5  (x, y)
extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
data source : C:\STATS\Lab_10\NDVI.nc 
names       : X2000.01.01, X2000.02.01, X2000.03.01, X2000.04.01, X2000.05.01, X2000.06.01, X2000.07.01, X2000.08.01, X2000.09.01, X2000.10.01, X2000.11.01, X2000.12.01, X2001.01.01, X2001.02.01, X2001.03.01, ... 
Date        : 2000-01-01, 2002-12-01 (min, max)
varname     : NDVI 

Question 3(c)

Aggregate all bands to generate a mean NDVI raster and a maximum NDVI raster; save the two new raster datasets to the hard drive.

mean(ndvi.rb)

class       : RasterLayer 
dimensions  : 360, 720, 259200  (nrow, ncol, ncell)
resolution  : 0.5, 0.5  (x, y)
extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
data source : in memory
names       : layer 
values      : 0, 0.8537528  (min, max)

writeRaster(mean(ndvi.rb),filename = "mean_ndvi.tif", overwrite=TRUE)
max(ndvi.rb)

class       : RasterLayer 
dimensions  : 360, 720, 259200  (nrow, ncol, ncell)
resolution  : 0.5, 0.5  (x, y)
extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
data source : in memory
names       : layer 
values      : 0, 0.9681  (min, max)

writeRaster(max(ndvi.rb),filename = "max_ndvi.tif", overwrite=TRUE)

Question 3(d)

Plot the maps of monthly NDVI of the year 2001

Question 3(e)

Create histograms of monthly NDVI of the year 2001.

hist(ndvi2001)

size changed to n because it cannot be larger than n when replace is FALSE39% of the raster cells were used (of which 42% were NA). 58093 values used.

Question 3(f)

Plot the NDVI map of July 2000; click any location with data on the map and retrieve the NDVI time series for all years; plot the NDVI time series of the selected location.

plot(ndvi.rb,7) 
values <- ndvi.rb[150, 50]
values <- click(ndvi.rb, n=1, xy=TRUE)

values <- click(ndvi.rb, n=1, xy=FALSE)